{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Draw Hopf fibration using python and POV-Ray"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What is Hopf fibration?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hopf fibration is a continous map from the 3-sphere $S^3$ onto the 2-sphere $S^2$, where the preimage of each point $p\\in S^2$ is a distinct circle called the **fiber** at $p$. The definition of the map is quite simple: identify $\\mathbb{R}^4$ with $\\mathbb{C}^2$ and $\\mathbb{R}^3$ with $\\mathbb{C}\\times\\mathbb{R}$ by writing $(x_1,x_2,x_3,x_4)$ as $(z_1,z_2)=(x_1+ix_2, x_3+ix_4)$ and $(x_1,x_2,x_3)$ as $(z,x)=(x_1+ix_2,x_3)$, thus $S^3$ is identified with the subset of $\\mathbb{C}^2$ such that $|z_1|^2+|z_2|^2=1$ and $S^2$ is identified with the subset of $\\mathbb{C}\\times\\mathbb{R}$ such that $|z|^2+x^2=1$, then the Hopf fibration is defined by\n",
    "$$(z_1,z_2)\\to (2z_1\\overline{z_2},\\, |z_1|^2-|z_2|^2).$$\n",
    "You can easily verify that if the left hand $(z_1,z_2)$ is on $S^3$, then the right hand must belong to $S^2$.\n",
    "\n",
    "It's also not hard to write down a parametric representation for the inverse map: a point $p=(x,y,z)\\in S^2$ can be parameterized by\n",
    "\n",
    "\\begin{align*}x &= \\sin(\\phi)\\cos(\\psi),\\\\ y &= \\sin(\\phi)\\sin(\\psi),\\\\ z &= \\cos(\\phi).\\end{align*}\n",
    "\n",
    "where $0 \\leq \\phi \\leq \\pi$ and $0 \\leq\\psi \\leq 2\\pi$. Then the fiber at $p$ is a circle on $S^3$ parameterized by $0\\leq\\theta\\leq2\\pi$:\n",
    "\n",
    "\\begin{align*}x_1&=\\cos((\\theta+\\psi) / 2)\\sin(\\phi / 2),\\\\ x_2&=\\sin((\\theta+\\psi) / 2)\\sin(\\phi / 2),\\\\x_3&=\\cos((\\theta-\\psi) / 2) \\cos(\\phi / 2),\\\\ x_4&=\\sin((\\theta-\\psi) / 2)\\cos(\\phi / 2).\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How can we visualize it?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To visualize the Hopf fibration we want to choose some points on the 2-sphere $S^2$, draw their fibers and and see what they look like. Since these fibers lie in the 4d space we cannot see them directly, but if we project them to 3d space using the [stereographic projection](https://en.wikipedia.org/wiki/Stereographic_projection) then some remarkable structure appears. The fibers are projected to circles in 3d space (one of which in a line, comes from the fiber through infinity), any two such circles are linked with each other and the line passes through all circles. The 3d space is filled with nested tori made of linking Villarceau circles, each tori is the preimage of a circle of latitude of the 2-sphere. \n",
    "\n",
    "So our plan is:\n",
    "\n",
    "1. Choose some points on the 2-sphere $S^2$.\n",
    "2. Compute their fibers as circles in $\\mathbb{R}^4$.\n",
    "3. Use stereographic projection to project these fibers to circles in 3d space and draw these circles.\n",
    "\n",
    "We will use [POV-Ray](www.povray.org) to render our 3d scene here. The computation task is handled in the python part and the rendering task is handled in the POV-Ray part. Certain background knowledge of POV-Ray's syntax is required to understand how the latter works. In summary we simply exports the data of the circles in the format of POV-Ray macros (\"macros\" are synonymous to \"functions\" in POV-Ray) and call these macros in the POV-Ray scene file. \n",
    "\n",
    "Some global settings:\n",
    "\n",
    "1. `POV_SCENE_FILE` is the scene file will be called by POV-Ray. It's named `hopf_fibration.pov` in the same directory with this notebook.\n",
    "2. `POV_DATA_FILE` is the output data file.\n",
    "3. `POV_EXE` is your POV-Ray executable file. (Don't forget add your POV-Ray executable file to system PATH!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import subprocess\n",
    "import numpy as np\n",
    "from IPython.display import Image\n",
    "\n",
    "PI = np.pi\n",
    "POV_SCENE_FILE = \"hopf_fibration.pov\"\n",
    "POV_DATA_FILE = \"torus-data.inc\"\n",
    "POV_EXE = \"povray\"\n",
    "COMMAND = \"{} +I{} +W500 +H500 +Q11 +A0.01 +R2\".format(POV_EXE, POV_SCENE_FILE)\n",
    "IMG = POV_SCENE_FILE[:-4] + \".png\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hopf inverse map and stereographic projection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hopf_inverse(phi, psi, theta):\n",
    "    \"\"\"Inverse map of Hopf fibration. It's a circle in 4d parameterized by theta.\n",
    "    \"\"\"\n",
    "    return np.array([np.cos((theta + psi) / 2) * np.sin(phi / 2),\n",
    "                     np.sin((theta + psi) / 2) * np.sin(phi / 2),\n",
    "                     np.cos((theta - psi) / 2) * np.cos(phi / 2),\n",
    "                     np.sin((theta - psi) / 2) * np.cos(phi / 2)])\n",
    "\n",
    "\n",
    "def stereo_projection(v):\n",
    "    \"\"\"Stereographic projection of a 4d vector with pole at (0, 0, 0, 1).\n",
    "    \"\"\"\n",
    "    v = normalize(v)\n",
    "    x, y, z, w = v\n",
    "    return np.array([x, y, z]) / (1 + 1e-8 - w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Circle passes through three points"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To draw the projected 3d circle of a fiber we choose three points on the fiber and construct the circle from their projected images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def normalize(v):\n",
    "    \"\"\"Normalize a vector.\n",
    "    \"\"\"\n",
    "    return np.array(v) / np.linalg.norm(v)\n",
    "\n",
    "\n",
    "def norm2(v):\n",
    "    \"\"\"Return squared Euclidean norm of a vector.\n",
    "    \"\"\"\n",
    "    return sum([x*x for x in v])\n",
    "\n",
    "\n",
    "def get_circle(A, B, C):\n",
    "    \"\"\"Compute the center, radius and normal of the circle passes\n",
    "       through 3 given points (A, B, C) in 3d space.\n",
    "       See \"https://en.wikipedia.org/wiki/Circumscribed_circle\"\n",
    "    \"\"\"\n",
    "    a = A - C\n",
    "    b = B - C\n",
    "    axb = np.cross(a, b)\n",
    "    center = C + np.cross((norm2(a) * b - norm2(b) * a), axb) / (2 * norm2(axb))\n",
    "    radius = np.sqrt(norm2(a) * norm2(b) * norm2(a - b) / (4 * norm2(axb)))\n",
    "    normal = normalize(axb)\n",
    "    return center, radius, normal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Convert vector/matrix to POV-Ray format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pov_vector(v):\n",
    "    \"\"\"Convert a vector to POV-Ray format.\n",
    "    \"\"\"\n",
    "    return \"<{}>\".format(\", \".join([str(x) for x in v]))\n",
    "\n",
    "\n",
    "def pov_matrix(M):\n",
    "    \"\"\"Convert a 3x3 matrix to a POV-Ray 3x3 array.\n",
    "    \"\"\"\n",
    "    return \"array[3]{{{}}}\\n\".format(\", \".join([pov_vector(v) for v in M]))\n",
    "\n",
    "\n",
    "# write a test to see if they work as expected:\n",
    "v = (1, 0, 0)\n",
    "print(\"POV-Ray format of {}: {}\".format(v, pov_vector(v)))\n",
    "M = np.eye(3)\n",
    "print(\"POV-Ray format of {}: {}\".format(M, pov_matrix(M)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Orient a circle in 3d space"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In POV-Ray CSG 3d circles are represented by the `Torus` object and by POV-Ray's default a `Torus` lies on the $xz$-plane with $y$-axis sticking through its center. So we need an orthogonal matrix to rotate it to a general orientation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def transform_matrix(v):\n",
    "    \"\"\"Return a 3x3 orthogonal matrix that transforms y-axis (0, 1, 0) to v.\n",
    "       This matrix is not uniquely determined, we simply choose one with a simple form.\n",
    "    \"\"\"\n",
    "    y = normalize(v)\n",
    "    a, b, c = y\n",
    "    if a == 0:\n",
    "        x = [1, 0, 0]\n",
    "    else:\n",
    "        x = normalize([-b, a, 0])\n",
    "    z = np.cross(x, y)\n",
    "    return np.array([x, y, z])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Export data to POV-Ray"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our POV-Ray macro as the interface between python and POV-Ray will be\n",
    "\n",
    "```\n",
    "Torus(center, radius, matrix, color)\n",
    "```\n",
    "\n",
    "This macro is implemented in the POV-Ray scene file. In python we just pack the data into this format and send them to POV-Ray for rendering."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def export_fiber(phi, psi, color):\n",
    "    \"\"\"Export the data of a fiber to POV-Ray format.\n",
    "    \"\"\"\n",
    "    A, B, C = [stereo_projection(hopf_inverse(phi, psi, theta))\n",
    "               for theta in (0, PI/2, PI)]\n",
    "    center, radius, normal = get_circle(A, B, C)\n",
    "    matrix = transform_matrix(normal)\n",
    "    return \"Torus({}, {}, {}, {})\\n\".format(pov_vector(center),\n",
    "                                            radius,\n",
    "                                            pov_matrix(matrix),\n",
    "                                            pov_vector(color))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Let's draw some examples!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we can draw a set of random points on $S^2$ and see how their fibers look like in 3d space: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_random_fibers(N):\n",
    "    \"\"\"Draw fibers of some random points on the 2-sphere.\n",
    "       `N` is the number of fibers.\n",
    "    \"\"\"\n",
    "    phi_range = (PI / 6, PI * 4 / 5)\n",
    "    psi_range = (0, 2 * PI)\n",
    "    phi_list = np.random.random(N) * (phi_range[1] - phi_range[0]) + phi_range[0]\n",
    "    psi_list = np.random.random(N) * (psi_range[1] - psi_range[0]) + psi_range[0]\n",
    "    with open(POV_DATA_FILE, \"w\") as f:\n",
    "        for phi, psi in zip(phi_list, psi_list):\n",
    "            color = np.random.random(3)\n",
    "            f.write(export_fiber(phi, psi, color))\n",
    "    subprocess.call(COMMAND, shell=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "draw_random_fibers(num_fibers=200)\n",
    "Image(IMG)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And also a flower pattern:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_flower(petals=7, fattness=0.5, amp=-PI/7, lat=PI/2, num_fibers=200):\n",
    "    \"\"\"parameters\n",
    "       ----------\n",
    "       petals: controls the number of petals.\n",
    "       fattness: controls the fattness of the petals.\n",
    "       amp: controls the amplitude of the polar angle range.\n",
    "       lat: controls latitude of the flower.\n",
    "    \"\"\"\n",
    "    with open(POV_DATA_FILE, \"w\") as f:\n",
    "        for t in np.linspace(0, 1, num_fibers):\n",
    "            phi = amp * np.sin(petals * 2 * PI * t) + lat\n",
    "            psi = PI * 2 * t + fattness * np.cos(petals * 2 * PI * t)\n",
    "            color = np.random.random(3)\n",
    "            f.write(export_fiber(phi, psi, color))\n",
    "    subprocess.call(COMMAND, shell=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "draw_flower()\n",
    "Image(IMG)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
